*cd "C:\Users\gianc\Downloads\Replication"
*global path "C:\Users\gianc\Downloads\Replication\"
cd "C:\Users\rodin\VenAnalytics Dropbox\VenAnalytics\Research\Papers\2025\Migration 2\Replication"
global path "C:\Users\rodin\VenAnalytics Dropbox\VenAnalytics\Research\Papers\2025\Migration 2\Replication"



*------------------------------------------
*Table 1: Monte Carlo Simulation Results Under No Cointegration
*------------------------------------------
clear all
set seed 5725

* Set number of replications and sample size
local reps = 1000
local nobs = 48

* Create matrix to store rejections
matrix results = J(`reps', 4, .) // cols: [1] levels, [2] diffs

forvalues i = 1/`reps' {
    
    clear
    qui set obs `nobs'

    * Create time variable
    qui gen t = _n
    qui tsset t

    * Generate independent white noise
    qui gen e1 = rnormal()
    qui gen e2 = rnormal()

    * Create random walks
    qui gen y = .
    qui gen x = .
    qui replace y = e1 in 1
    qui replace x = e2 in 1

    forvalues j = 2/`nobs' {
        qui replace y = y[_n-1] + e1 in `j'
        qui replace x = x[_n-1] + e2 in `j'
    }

    * Engle-Granger on levels
    quietly egranger y x
	local Z_levels=e(Zt)
    local p_levels=0
	if e(Zt)<e(cv5) {
		local p_levels = 1
		} 
	

    * First differences
    qui gen dy = D.y
    qui gen dx = D.x
    quietly egranger dy dx
	local Z_diffs=e(Zt)
    local p_diffs=0
	if e(Zt)<e(cv5) {
		local p_diffs = 1
		} 
	
	display `i'

    * Store results
    matrix results[`i',1] = `Z_levels'
    matrix results[`i',2] = `p_levels'
    matrix results[`i',3] = `Z_diffs'
    matrix results[`i',4] = `p_diffs'
}

* Calculate mean rejection rates
svmat results

*Table 1
tabstat results*, stats(mean sd min max) save

*------------------------------------------
* Figure 1: Distribution of Engle-Granger Test Statistics Using Level Versus First-Differenced Data
*------------------------------------------
 twoway (histogram results1, percent color(gs12%50) lwidth(none)) ///
       (histogram results3, percent color(blue%50) lwidth(none)), ///
       xline(-3.472, lcolor(black%100) lpattern(dash) lalign(inside)) /// 	   
	   text(0.5 -3.3 "5% critical value", place(e) size(small)) ///
legend(order(2 "Bahar-Hausmann (First Differences)" 1 "Engle-Granger (Levels)" ) rows(1) position(6) size(small)) ///
	  	   text(13.5 -6 "Test finds evidence of cointegration", size(small)) ///
       text(13.5 -1.2 "Test finds no evidence of cointegration", size(small)) ///
	    ylabel(, nogrid) xlabel(, nogrid) xscale(range(-11 1.5)) saving("Figure/hist.gph", replace)  plotregion(margin(zero))
 graph export "Figure/hist.jpg", replace
*------------------------------------------
 
*------------------------------------------
*Table 2 and Appendix
*------------------------------------------
*Data preparation
*------------------------------------------
  

import delimited "${path}/BH_replication_package/data/DCOILBRENTEU.csv", delimiter(",") clear
destring dcoilbrenteu , replace force
collapse (mean) brentprice=dcoilbrenteu, by(month year)
tempfile price
save `price'

import delimited "${path}/BH_replication_package/data/venezuela_crude_MOMR_OPEC.csv", delimiter(",") clear
ren production oil_productionk
g oil_production1m = oil_productionk/1000
tempfile oil
save `oil'

import delimited "${path}/BH_replication_package/data/Data_NSA.csv", delimiter(",") clear
drop oil_production*
drop if mi(month) | mi(year)
merge 1:1 month year using `oil', keep(master matched) nogen
merge 1:1 month year using `price', keep(master matched) nogen
g oil_incomeM = 30*oil_production1m*brentprice
drop if mi(swencounters_vzla_nsa)


g date_str = string(month)+"/"+string(year)
gen date=monthly(date_str, "MY")
format %tm date

tsset date, monthly

tempfile g1 g2 g3 g4

*log variables
g log_encounters_nsa = log(swencounters_vzla_nsa)
g log_income = log(oil_incomeM)
g log_price = log(brentprice)
g log_production = log(oil_productionk)

*Results matrix
mat A=J(21,3,.)
mat B=J(21,3,.)
mat C=J(15,4,.)
mat D=J(15,4,.)
*------------------------------------------

*Bahar and Hausmann (2025a, 2025b) test statistic
egranger d.swencounters_vzla_nsa d.oil_incomeM

*Table 2: Engle-Granger Tests for Cointegration, Logarithmic Transform
*------------------------------------------
tokenize log_income log_price log_production

forvalues i=1/3  {
*Engle-Granger
    egranger log_encounters ``i''
    mat A[1,`i'] = e(Zt)
    mat A[2,`i'] = e(N)
    mat A[3, `i']=e(cv1)
    mat A[4, `i']=e(cv5)
    mat A[5, `i']=e(cv10)
	predict res, resid
	qui wntestb res
	mat A[6, `i']=r(stat)
	mat A[7, `i']=r(p)
	drop res
*Aug. Engle-Granger (12 lags)
 egranger log_encounters ``i'', lags(12)
    mat A[8,`i'] = e(Zt)
    mat A[9,`i'] = e(N)
    mat A[10, `i']=e(cv1)
    mat A[11, `i']=e(cv5)
    mat A[12, `i']=e(cv10)
	predict res, resid
	qui wntestb res
	mat A[13, `i']=r(stat)
	mat A[14, `i']=r(p)
	drop res
*Aug. Engle-Granger (12 lags + trend) 
egranger log_encounters ``i'', lags(12) trend	
    mat A[15,`i'] = e(Zt)
    mat A[16,`i'] = e(N)
    mat A[17, `i']=e(cv1)
    mat A[18, `i']=e(cv5)
    mat A[19, `i']=e(cv10)
	predict res, resid
	qui wntestb res
	mat A[20, `i']=r(stat)
	mat A[21, `i']=r(p)
	drop res
}
*------------------------------------------

*Table A.1: Bahar-Hausmann Implementation of Engle-Granger Test
*------------------------------------------
egranger d.swencounters_vzla_nsa d.oil_incomeM, reg
*------------------------------------------

*Table A.2: Bahar-Hausmann Implementation of Engle-Granger Test
*------------------------------------------
tokenize oil_incomeM brentprice oil_production1m

forvalues i=1/3  {
*Engle-Granger
    egranger swencounters_vzla_nsa  ``i''
    mat B[1,`i'] = e(Zt)
    mat B[2,`i'] = e(N)
    mat B[3, `i']=e(cv1)
    mat B[4, `i']=e(cv5)
    mat B[5, `i']=e(cv10)
	predict res, resid
	qui wntestb res
	mat B[6, `i']=r(stat)
	mat B[7, `i']=r(p)
	drop res
*Aug. Engle-Granger (12 lags)
 egranger swencounters_vzla_nsa  ``i'', lags(12)
    mat B[8,`i'] = e(Zt)
    mat B[9,`i'] = e(N)
    mat B[10, `i']=e(cv1)
    mat B[11, `i']=e(cv5)
    mat B[12, `i']=e(cv10)
	predict res, resid
	qui wntestb res
	mat B[13, `i']=r(stat)
	mat B[14, `i']=r(p)
	drop res
*Aug. Engle-Granger (12 lags + trend) 
egranger swencounters_vzla_nsa  ``i'', lags(12) trend	
    mat B[15,`i'] = e(Zt)
    mat B[16,`i'] = e(N)
    mat B[17, `i']=e(cv1)
    mat B[18, `i']=e(cv5)
    mat B[19, `i']=e(cv10)
	predict res, resid
	qui wntestb res
	mat B[20, `i']=r(stat)
	mat B[21, `i']=r(p)
	drop res
}
*------------------------------------------


*Table A.3: Dickey-Fuller and Augmented Dickey-Fuller Tests, Logarithmic Transform
*------------------------------------------
tokenize log_encounters log_income log_price log_production

forvalues i=1/4  {
    dfuller ``i''
    mat C[1,`i'] = r(Zt)
    mat C[2,`i'] = r(N)
    mat C[3, `i']=r(cv_1)
    mat C[4, `i']=r(cv_5)
    mat C[5, `i']=r(cv_10)
 
 dfuller ``i'', lags(12)
    mat C[6,`i'] = r(Zt)
    mat C[7,`i'] = r(N)
    mat C[8, `i']=r(cv_1)
    mat C[9, `i']=r(cv_5)
    mat C[10, `i']=r(cv_10)

dfuller ``i'', lags(12) trend	
    mat C[11,`i'] = r(Zt)
    mat C[12,`i'] = r(N)
    mat C[13, `i']=r(cv_1)
    mat C[14, `i']=r(cv_5)
    mat C[15, `i']=r(cv_10)
}
*------------------------------------------


*Table A.4: Dickey-Fuller and Augmented Dickey-Fuller Tests, Untransformed Series
*------------------------------------------
tokenize swencounters_vzla_nsa oil_incomeM brentprice oil_production1m

forvalues i=1/4  {
    dfuller ``i''
    mat D[1,`i'] = r(Zt)
    mat D[2,`i'] = r(N)
    mat D[3, `i']=r(cv_1)
    mat D[4, `i']=r(cv_5)
    mat D[5, `i']=r(cv_10)
 
 dfuller ``i'', lags(12)
    mat D[6,`i'] = r(Zt)
    mat D[7,`i'] = r(N)
    mat D[8, `i']=r(cv_1)
    mat D[9, `i']=r(cv_5)
    mat D[10, `i']=r(cv_10)

dfuller ``i'', lags(12) trend	
    mat D[11,`i'] = r(Zt)
    mat D[12,`i'] = r(N)
    mat D[13, `i']=r(cv_1)
    mat D[14, `i']=r(cv_5)
    mat D[15, `i']=r(cv_10)
}
*------------------------------------------
*------------------------------------------


*------------------------------------------
*Summary of results
*------------------------------------------
*Table 2
svmat A 
*Table A.2
svmat B 
*Table A.3
svmat C 
*Table A.4
svmat D

keep A* B* C* D* 
save "${path}/Results/results.dta"
*------------------------------------------





